Structure of molecular liquids: cavity and bridge functions of the hard spheroid fluid 
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We present methodologies for calculating the direct correlation function, c(l,2), the cavity func- 
tion, 2/(1,2), and the bridge function, fe(l,2), for molecular liquids, from Monte Carlo simulations. 
As an example we present results for the isotropic hard spheroid fluid with elongation e = 3. The 
simulation data are compared with the results from integral equation theory. In particular, we solve 
the Percus-Yevick and Hypernetted Chain equations. In addition, we calculate the first two terms 
in the virial expansion of the bridge function and incorporate this into the closure. At low densi- 
ties, the bridge functions calculated by theory and from simulation are in good agreement, lending 
support to the correctness of our numerical procedures. At higher densities, the hypernetted chain 
results are brought into closer agreement with simulation by incorporating the approximate bridge 
function, but significant discrepancies remain. 
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PACS numbers: 61.20.Gy, 61.20.Ja, 05.20.Jj 



I. INTRODUCTION 



The equilibrium properties of homogeneous fluids of 
spherical particles have been extensively studied both by 
theory and simulation and a great deal is now known 
about the thermodynamic properties and the fluid struc- 
ture P, Q ■ Simulation has been used to calculate the 
total and direct correlation functions |3| ithe cavity func- 
tion 0, and the bridge function 0, 13 13 • On the theo- 
retical side, integral equation theory (IET) is now capable 
of making some very accurate predictions. Percus-Yevick 
(PY) and Hypernetted Chain (HNC) theories have now 
been extended, for example, by mixing closures so as to 
obtain identical virial and compressibility equations of 
state 0, 0] . An alternative approach has been to incor- 
porate approximate forms for the bridge function in the 
HNC closure. These may take the form of a low-order 
virial expansion, a bridge function from a reference fluid 
or an approximate closure relation 0, ll2L Il3| . While 
there is still work to be done, especially perhaps on a 
fundamental treatment of the bridge function, the foun- 
dations are rather solid. As a consequence, one is in a 
good position to construct good density functionals to 
describe inhomogeneous fluids. A key ingredient of a 
density functional is an assumed form for the inhomoge- 
neous direct correlation function. It is clearly reassuring 
if this quantity reduces to the known homogeneous func- 
tion in the uniform limit and, for spherically symmetric 
particles, this is a test one may apply. 

The equilibrium properties of isotropic fluids of non- 
spherical particles are less well-characterised. There are 
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relatively few simulation studies on the direct correlation 
function 0] and little data for the cavity or bridge func- 
tions have been published (site-site functions have been 
calculated for hard sphere dirners and water 0, 0, 0] 
while the first bridge diagram has been calculated for 
the hard spherocylinder fluid for a number of fixed orien- 
tations The PY and HNC equations have been 
solved for axially symmetric particles (e.g. hard ellip- 
soids, hard spherocylinders and truncated hard spheres) 
and the general conclusion is that HNC is superior to 
PY for significantly aspherical particles, but that there 
is still a substantial discrepancy between theory and sim- 
ulation, especially at high density |l9U2CU2lt l22 | . There 
have been some attempts to go beyond HNC. Pospisil 
et al. [2j| have investigated the use of a modified Verlet- 
bridge closure and have reported improved results. Singh 
et al. 0| applied a non-spherical version of the Rogers- 
Young method of mixing PY and HNC closures, again 
obtaining results for spheroids in good agreement with 
simulation. Nevertheless the number of such studies is 
relatively small and, as yet, we do not have sufficient 
simulation and theoretical studies to claim a foundation 
to rival that enjoyed by spherical particles. 

In this paper we try to address some of these issues, 
using Monte Carlo (MC) simulations and IET. On the 
simulation front, we present methodologies for calculat- 
ing the direct correlation, cavity and bridge functions 
for isotropic fluids of axially symmetric particles using 
advanced MC techniques. These methods are used to 
calculate the molecular correlation functions for a fluid 
of hard spheroids with major axis of length a and mi- 
nor axis of length b. We focus here on an elongation 
e = a I b = 3, and present results for a range of den- 
sities in the isotropic phase. IET is adapted for fluids 
of anisotropic particles using invariant expansions of the 
correlation function s 1251 1261 a nd efficient numerical al- 
gorithms [2(], I27], l2St l29t l30l . In particular we use the 
relaxation method of Ng [3l| to provide a robust and 



2 



easily-programmable algorithm for numerically solving 
the integral equations. We also examine some analytical 
properties of the cavity function for non-spherical hard 
particles and calculate the first two terms in the virial 
expansion of the bridge function. 

The paper is organised as follows. In section [H] we 
give the basic equations relating the correlation functions 
studied in this article. Section ITTT1 describes the simula- 
tion methods used for the calculation of the cavity func- 
tion and bridge function. In Section ll VI we present some 
technical details of the numerical solution of the IET us- 
ing the method of Ng j3]| and the Monte Carlo proce- 
dure used to compute the bridge diagrams. In Section 
El the results of simulations and IET are compared and 
discussed. The general conclusion of our study and some 
future avenues of work are given in Section IVII Some 
more technical details of the Monte Carlo algorithm for 
the calculation of the cavity function are presented in the 
Appendix. 



II. GENERAL FORMALISM 

The structure of a fluid may be described at a two- 
particle level by the total correlation function (TCF) 
h(l,2) = g(l,2) — 1 (where .9(1,2) is the pair distribu- 
tion function) or the direct correlation function (DCF) 
c(l,2). These are linked via the Ornstein-Zernike (OZ) 
equation, which for a homogeneous fluid of axially sym- 
metric molecules is 0,0 

6.(1,2) = c(l,2) + |- J d3c(l,3)M3,2) (1) 

where p is the number density and, as is traditional, (i) — > 
(ri,Ui). Here denotes the centre of mass position of 
particle i whilst itj represents a unit vector along the 
particle's symmetry axis. 

To determine h(l, 2) and c(l, 2), Eq. Q is usually sup- 
plemented by an approximate closure relation. These 
take the form 

c(l,2) = (l + ft(l,2))(l-exp[^(l,2)]) PY (2) 

c(l, 2) = h(l, 2) - log[l + h{l, 2)] - /3V(1, 2) HNC 

(3) 

where V(l,2) is the intermolecular potential and (3 = 
l/k B T. 

The exact closure relation can be written as follows [3 



1/(1, 2) = cxp {h(l, 2) - c(l, 2) + 6(1, 2)} (4) 

where 6(1, 2) is the bridge function and y(l, 2) is the cav- 
ity or background correlation function defined by the re- 
lation: 

1/(1,2) = 5 (l,2)exp[/3V(l > 2)] . (5) 



Eq. may be regarded as a definition of 6(1, 2) and the 
approximate closure relations may be regarded as ap- 
proximations to the unknown 6(1,2). In particular, the 
PY and HNC closures, Eqs. 10,© respectively, corre- 
spond to 

6(l,2) = 7 7 (l,2)-log(l + 7 7 (l,2)) PY (6a) 
6(1, 2) = HNC (6b) 

where 77(1, 2) = h(l, 2) - c(l, 2). 

The bridge function may be expressed as a virial ex- 
pansion 

6(1,2) = (1,2) , 

n>2 

where 5B„(1, 2) are the bridge diagrams. In principle, this 
provides a route for the exact calculation of the bridge 
function, but in practice it is only feasible to calculate 
low-order terms, as has been done for hard spheres 0, 
Il3j . In this paper we use the two lowest-order estimates 
of the bridge function for hard spheroids 

62(1,2) =p 2 Q5 2 (l, 2) HNC+B2 (6c) 

63(1,2) = (0 2 <B 2 (l,2) + p 3 » 3 (l,2) HNC+B3 (6d) 

to extend the HNC closure relation. In this paper we 
investigate all four closure relations, i.e. Eq. Q with the 
bridge function specified by one of Eqs. (|5a}-|(63J|. 

The numerical solution of the integral equation and 
MC calculation are based upon the expansion of two- 
particle functions in a basis set of rotational invariants 

urn 

F(l,2) = Y i F mnt {r)^ mnt {u 1 ,u 2 ,u r ) , (7) 

rant 

*"^(« 1 ,ti 2> « r ) = 47r y ( m n M 

^ ' ' ^ ^ \Xl X2 Xr 

x Y mxl (u 1 )Y nX2 (u 2 )Ce Xr (u r ) , (8) 

where r is the intermolecular distance; u r is a unit vector 
along the intermolecular vector, U\, u 2 are the orienta- 
tions of the molecules in a given system of coordinates 
('laboratory frame'), Y mx (u) are the spherical harmonics 
functions, C mx (u) = (An/(2m+ l)) 1 / 2 Y,„ x (u) and 

I m n £ \ 

V Xl X2 Xr J 

are the standard 3j symbols. 

Some quantities of interest are easier to compute in a 
system of coordinates whose z-axis lies along the inter- 
molecular vector ('molecular frame'). The expansion in 
the molecular frame has the form: 

F(l, 2) = 47T f m nxW^("l)^(«2) , (9) 
mnx 
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where x = ~X- The two sets of coefficients are connected 
through the ^-transform and its inverse: 



Fmnx( r ) — ^ ] 



X X o 



F mn \r) 



(10a) 



F™^(r) = (2£ + l)J2{™ I J ) /7 ""^ ,ri ' ,1 " 1>1 



III. SIMULATION METHOD 

A. Direct correlation function 

The total correlation function may be determined di- 
rectly from simulation through the pair distribution func- 
tion g(l,2). The spherical harmonic coefficients are de- 
termined as usual from 

111 



9mn X (r) = 47T ff0 ooW (Y^ x (ui)Y*Ju 2 )) 7 



(11) 



where Y mx {u) is a spherical harmonic, x — ~Xi 9ooo( r ) is 
the pair distribution function of the particle centres, and 
the angled brackets denote an average over all molecules 
in the shell [r, r + Sr}. These coefficients are defined in 
the molecular frame described in Section ^ |2j ■ From 
Eq. JTTJ it follows that h mnx {r) = g mnx (r) - S mQ 6 n0 6 x o. 

The direct correlation function may be found from the 
measured total correlation function in two ways. In re- 
ciprocal space, using the molecular frame expansion, the 
Ornstein-Zcrnike equation becomes 



'rnnx 



(*0 



x (k) + (-l)Xp V h mjx (k)c jnx (k) (12) 



j 



where f(k) is the (three-dimensional) Fourier transform 
of a function f(r). The structure of this equation with 
respect to the first two indices leads to a matrix notation 



HJk) = C x (k) + (-ir P HJk)CJk) 



(13) 



where Pf{x) — x~ 1 dPi(x)/dx and Pi(x) are Legendre 
polynomials. A x-transform, Eq. (|10a|) . then converts 
the functions to the required molecular frame. It is as- 
sumed that a separation R exists such that Q x (r) = and 
C x (r) = for all r > R. Equation i|14bjl is solved itera- 
tively to find Q x (r) from the functions H x (r) determined 
in the simulation. Once this procedure has converged, 
Eq. (|14a|) is used to determine C x (r). At very small r 
this involves the difference between two large quantities 
possibly leading to numerical difficulties. These may be 
avoided by a procedure outlined previously finding 
C x (0) from 



C x (0) - fL(0) = (-l)*27rp 




dr rH x (r)Q' x 



r) 



QUr)Qj'(r))-Q x (0)QT(0) (16) 



c(l, 2) may also be obtained via a real-space factoriza- 
tion [3, |33, |34j . It is possible to write 



rt x (r) = -Q' x (r) + 2tt(-1)*p / dsQ' x (s)Ql(s - r) 

J r 

(14a) 

f R 

rH x (r) = -Q'Jr) + 27r(-l) x p / ds(r - s)H' (r - s)Q x {s) 

Jo 

(14b) 

where the new matrix Q x (r) has been introduced, 
Q' x (r) — dQ x (r)/dr, and Q x (r) is the transpose of Q x (r). 

The so-called 'hat' transform giving the functions H x , 
C x , that appear in Eqs. (fHjl is defined in the laboratory 
frame 



mni j 



ds s- l r nt {s)P!{r/s) (15) 



and the C x (r) for r — > are then determined by interpo- 
lation. 



B. Cavity correlation function 

There are two methods for the calculation of the cavity 
function, Eq. JSJ, either by a direct simulation of two non- 
interacting cavity particles Q| or through the test-particle 
method based on Henderson's equation |35| . The first of 
these methods is more useful for large cavity separations 
while the second is better as r — > 0. 



1. Direct simulation method 

The direct simulation method follows from the obser- 
vation that for the hard particle fluid, the cavity function 
may be identified as the pair distribution function for a 
pair of non-interacting cavities [3(|. In a MC simula- 
tion it is convenient to constrain the two cavities to be 
within a given range of separations r*i2 • Even so, in a nor- 
mal MC simulation the probability distribution P C av(l, 2) 
is likely to vary rapidly with separation r^, leading to 
poor sampling in the regions where the function is rela- 
tively small. To circumvent this problem, the umbrella 
sampling technique is employed |37|. The r-separation 
of the cavities is divided into a set of overlapping win- 
dows. Within each window a weight function w{rvi) is 
introduced into the Monte Carlo moves; this function is 
iteratively refined so as to produce a flat sampled proba- 
bility distribution. This weight may be subsequently re- 
moved to give the true probability distribution for each 
window and the full distribution is recovered using the 
self-consistent histogram method [38| . 

The cavity function is, to within a multiplicative con- 
stant, equal to P cav (l,2)/rf 2 . When the cavity particles 
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are constrained this constant cannot be determined di- 
rectly . However it may be found by enforcing the con- 
dition that y(l, 2) = 1 + h(l, 2) when outside the overlap 
region 0|. 

For this scheme to be effective a good choice of the 
weighting function is needed. For hard spheres a good 
choice proved to be an analytic approximation to y(r) 
Here we employ a more general method based on the 
Wang-Landau method [Sfj, Efij . Briefly, this is an itera- 
tive method that updates an initial guess to the weight 
function using a decreasing modification factor. Full de- 
tails are given in the appendix. The implementation used 
here is similar in spirit to the extended density of states 
method (EDOS) The spherical harmonics coef- 

ficients y m nx( r ) are found in the same way as those for 
the pair distribution function g m nx( r )- 



by 



The spherical harmonic coefficients y m nx( r ) are given 



ymn X ( r ) = (4tt) 1 / duodux y(0,l)Y*(u o )Y*Ju 1 ) 



N 



= (4 7 r)- 1 exp[/3 Mcx ] ^exp I -^/?y(0,j) 

Xy m X N^(«l)) , 



(18) 



where the angled brackets denote averages over test par- 
ticle insertions in the range [r, r + Sr]. 



C. Bridge function 



2. Test particle method 

In the canonical ensemble, Henderson's equation for a 
system containing N molecules may be written [f| 



y(0, 1) = exp [flue 



N 



exp -J2/3(V(0,j) 



(17) 



J>2 



N,V.T 



where /i cx is the excess chemical potential and the an- 
gled brackets denote an ensemble average over particles 
1, . . . , N. The term in the angled brackets corresponds to 
the Boltzmann factor of a molecule with the interaction 
with another molecule 1 neglected. This may intuitively 
be equated to a fluid consisting of 2 non-interacting cav- 
ities and N — 1 other molecules. Additionally this is also 
equivalent to the calculation of the acceptance criterion 
in a Metropolis MC simulation - the interaction between 
a hypothetical molecule with position t*o and orientation 
Uq and every molecule in the system apart from 1 is the 
quantity that is calculated when an attempt is made to 
move molecule 1 to position rg and orientation Uq. So as 
any molecule in the system may be labelled 1 the quan- 
tity in the angled brackets in Eq. l(T?|) is calculated for 
every attempted MC move. This fact has been used in 
previous studies of atomic fluids |y, |7fl allowing the calcu- 
lation of 2/(0, 1) at essentially no extra cost. However, for 
molecular fluids, where the maximum angular displace- 
ment in a Monte Carlo simulation may be much smaller 
than 2n, this would lead to poor sampling of the angular 
dependence y(Q, 1). Hence the calculation of y(Q, 1) pro- 
ceeds by inserting a number of test particles (labelled 0) 
in the vicinity of each molecule in the simulation (labelled 
1). The Boltzmann factor (neglecting the interaction be- 
tween and 1) is then calculated. For hard molecules, 
as in the present case, this is simply if the test particle 
overlaps with any other molecule (excluding molecule 1) 
or 1 if there are no overlaps. 



Once spherical harmonic expansions for /i(l,2), c(l,2) 
and 2/(1,2) have been determined, the final step is to 
invert Eq. (@J for 6(1, 2). While the presence of the expo- 
nential on the right hand side of Eq. {3J is troublesome 
for the spherical harmonics expansions, it may be easily 
circumvented |19|. Taking the logarithm and differenti- 
ating Eq. Q| with respect to r gives 



dy(l,2) 
dr 



2/(1,2) 



dh{l,2) <9c(l,2) 96(1,2) 



Or 



dr 



dr 



Inserting the spherical harmonic expansions of the pair 
functions and integrating over angles gives 

dy m n X ( r ) _ a r mm'm" r nn'n" ( \ 

d r — m Z^i XX'x" 1 xx'x"» m ' " X"V ) 



m n x 
m n" x 



(20) 



where 



T^x™ = duY^ x {u)Y m , x ,{u)Y m „ x „{u) 



'(2to' + 1)(2to" + 1) „ , 



y 47r(2m + l) 
x C(m",m',m;x",x',X) » 



C(m",m',m; 0,0,0) 
(21) 



and where C(m", to', to; x" \ x'x) are Clebsch-Gordan co- 
efficients. Eq. (I20[) can be solved using standard numer- 
ical methods |43| for the derivatives db mnx (r)/dr, and 
these are integrated numerically to give the bridge func- 
tion components b mnx (r). 



D. Simulated system 

The simulated system consists of a fluid of hard pro- 
late spheroids of elongation e = a/6 = 3. This is a com- 
mon model for molecular fluids and liquid crystals and 
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along with similar models such as hard spherocylinders 
has been well studied |44|. 

For the calculation of h(l,2), systems of 2048 
molecules were simulated using constant NVT MC sim- 
ulations. Data for the calculation of h(l,2) were gath- 
ered every 500 MC sweeps (each sweep is on average 
1 attempted translation and 1 attempted rotation per 
molecule) over a total of 5 x 10 5 MC sweeps. The c mnx (r) 
coefficients were then calculated from the h mnx (r) co- 
efficients following Sec. 1111 Al The spherical harmon- 
ics expansions for the pair functions were truncated at 
wi max = n max = 8 and the grid spacing Sr = 0.016. 

For the calculation of y(l, 2), systems of 512 molecules, 
including 2 cavity molecules, were simulated (smaller sys- 
tems are sufficient for the calculation of y(l,2) as its 
long-range behaviour is identical to h{\, 2)). The r sep- 
aration between the cavity particles was split into over- 
lapping windows covering r/b = [0.03,0.50], [0.20,1.20], 
[1.00,2.00], [1.80,2.80], and [2.60,3.60]. In each window 
the weight function was determined over at least 15 iter- 
ations (see appendix for details). Once the final weight 
function was determined, 7/(1, 2) data were gathered over 
a total of 2 x 10 7 MC sweeps. Error estimates were 
made by splitting this into 4 subruns. 2/(1,2) was cal- 
culated for the region r/b — [0,0.15] using the test par- 
ticle insertion method (Sec. IIIIBH) . 6.(1,2) and y(l,2) 
have been calculated at reduced densities p* — p/ p cp — 
0.10,0.20,0.30,0.40,0.50 where p cp = V2/(ab 2 ) is the 
close-packed density. 



IV. INTEGRAL EQUATION THEORY 

To solve the integral equations in the isotropic phase 
we have used the standard rotationally invariant decom- 
position of the angular part of the correlation functions 
as discussed in detail in Refs. [2(1 |2lJ • The solution is 
calculated iteratively with the help of the method of Ng 
|3l| that yields fast convergence even at densities close 
to those where no real solu tion exists. 

We describe in short the iNd method as ap plie d to the 
hard spheroid fluid. An iteration step in the INd method 
is done using as input a linear combination of the p func- 
tions obtained in the p previous steps. The coefficients 
of the linear combination are calculated from a smallest 
displacement condition. 

An iteration has the generic form 



/ i+1 (l, 2) =0^(1,2)] 



(22) 



with ti(l, 2) = Ml, 2)-J2 ai, m A/ i , m (l, 2) (23) 



where 0[f] is the iteration operator, /,(1,2) is the 
iteration result and A/j, m (l,2) = /i(l,2) - /i_ m (l,2). 
At each iteration step the scalars a,i >m are computed from 
the minimum condition of the following functional: 



Close to the solution we assume that the differences 
A/i,ro(l>2) are small and we expand Eq. ljL'21) up to the 
first order: 



0[/i(l,2)- ai, m Afi, m (l,2)] 

771—1 

» 0[/i(l,2)] - J2 V ^fflf ^tM) • (25) 

The coefficients «i. m that satisfy the approximated min- 
imum condition, Eq. (|24H . are the solutions of a linear 
system of equations 



^ akmOLi^n = h , k = 1. ..p , 



(26) 



m— 1 



where the coefficients a km and bk are determined from 
the following equations: 



d2[f l+1 (l,2)-t l (l,2)Y 



(24) 



a km = J d2 SO[fi(l, 2)} k SO[fi(l, 2)] m (27) 
b k = J d2 (0[/i(l,2)] - fi(l,2))SO[fi(l,2)] k , (28) 

and 

50[/41,2)] fc = ^M^l A/ . fc(1)2) . 
In our case the nonlinear operator 0[-} has the form 

m 

0[rj\ = A(l, 2)(-l - 77(1, 2)) + (1 - A(l, 2))c cl (l, 2) , (29) 

where A(l, 2) has the value 1 if spheroids 1 and 2 overlap 
and the value if they do not. c c i(l, 2) is given by 

PY (30a) 

exp(ry(l, 2)) - 77(1,2) - 1 HNC (30b) 

exp(?7(l, 2) +62(1,2)) -77(1,2) - 1 HNC+B2 (30c) 

exp(77(l,2) + 6 3 (l,2)) -77(1,2)- 1 HNC+B3 (30d) 

corresponding to the closure relations of Eqs. 0. 

We mention that the indirect correlation function, 
77(1, 2) = 6,(1, 2) - c(l, 2), that appears in Eqs. ^ 
is computed at each iteration step from the OZ equation, 
and the expansion l|25|) is performed with the function 
?7(1,2). 

The algorithm is written using the angular components 
of the operator 0[), Eq. I|29|l . and of the correlation func- 
tions c(l, 2), 77(1, 2) as described in full detail in Ref. [20| . 

In the numerical calculation, the expansion in rota- 
tional invariants of the correlation functions, Eq. Q, is 
truncated at m max = ?i max — 8 and all non-zero com- 
ponents consistent with this truncation are kept. The 
integral equation was discretized on a grid in steps of 
0.016. 
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The first- and second-order bridge diagrams were com- 
puted using an extension of the Monte Carlo methods 
described by Ree and Hoover 0, lift The first step was 
to convert the diagrams from Ref . [l2j > given in terms of 
Mayer /-bonds, into Ree-Hoover diagrams, where field 
points are connected either by an /-bond or by an e-bond, 
where e = 1 + /. The overall bridge function is obtained 
from a weighted sum of these Ree-Hoover diagrams. Par- 
ticle 1 is placed at the origin with its symmetry axis along 
the z axis, and a second particle is placed at random so 
that it overlaps the first particle. A third particle is sim- 
ilarly randomly placed to overlap the second particle and 
so on. When calculating the first set of bridge diagrams, 
a chain of four such particles is generated. The second set 
of bridge diagrams require a five-particle chain. The over- 
laps between all pairs of particles are checked. If the con- 
figuration corresponds to one of the Ree-Hoover bridge 
diagrams, then the separation between the two end par- 
ticles of the chain is calculated, ready for accumulation 
as a histogram. To obtain the angular expansion coef- 
ficients, the Ree-Hoover weighting is multiplied by the 
spherical harmonic product Y*(ui)Y*^(u2) , where the 
unit vectors are expressed relative to the vector joining 
the two end particles of the chain (1 and 2). The com- 
ponents of the bridge diagrams are thus in the molecular 
frame. After a sufficient number of Monte Carlo config- 
urations have been generated, N con {, the final results for 
the bridge function are obtained by normalising the his- 
togram values, firstly by a factor of iV con f, secondly by a 
factor of the volume of the spherical shell corresponding 
to the separation between the particles and thirdly by 
an appropriate power of the pair excluded volume (i.e. 
the square for the first-order term and the cube for the 
second-order term) . Errors may be estimated in the stan- 
dard way, by dividing the total number of configurations 
into sub-batches and calculating sub-averages. 

We used 1.6 x 10 9 trial chain configurations to obtain 
the first bridge diagram and 1.1 x 10 9 trial chain configu- 
rations to obtain the second bridge diagram. The relative 
error estimate is close to 1% except for the r < 0.16 do- 
main that is sampled poorly by this method. 

In summary we now have four sets of integral equation 
results with which to compare simulation, correspond- 
ing to the four closures of Eqs. ©, namely PY, HNC, 
HNC+B2 (first-order bridge) and HNC+B3 (first- and 
second-order bridges). 



V. RESULTS 

A. Equation of state and stability with the bridge 
diagrams 

The angular coefficients of the direct and the indirect 
correlation functions obtained from PY and HNC inte- 
gral equations for non-spherical particles have already 
been extensively compared with simulation results in 
Refs. HHmEl. We limit our discussion to the effect 
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FIG. 1: The direct correlation function components 000 and 
220 (in the lab frame) at reduced densities p* =0.1 (left) and 
p" = 0.5 (right) obtained from IET and MC. 
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FIG. 2: The equation of state (upper panel) obtained by the 
virial (v) and compressibility (c) routes; and the Kerr coef- 
ficients (lower panel) obtained from pure HNC equation and 
HNC with the first two bridge diagram corrections. The lines 
for the Kerr coefficients correspond to m = 2, 4, 6, 8 in ascend- 
ing order and the symbols represent the MC data (equation 
of state data from Ref. H^].) 

of the inclusion of the bridge diagrams in the closure. 
Fig. for two angular components of the direct correla- 
tion function at two densities, shows that the agreement 
between MC data and IET improves at high density if 
the HNC closure is supplemented by the inclusion of the 
low-order bridge diagrams. 

We have a mixed picture for the equation of state and 
Kerr or stability coefficients. The latter gives a measure 
the stability of the isotropic phase relative to the nematic 
phase. The isotropic phase is stable when 0> US 

1 - (2m + l)- 1 / 2 ?" 1 ™ ^) > , m = 2,4,6,... (31) 

where c mm0 (0) is the low-fc limit of the Fourier- 
transformed direct correlation function component 
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FIG. 4: The same functions as in Fig.[3]at density p* — 0.3. 



c mm0 (r) in the laboratory frame. Fig. [2] shows that the 
inclusion of the first-order bridge diagram improves the 
agreement with the MC data for both the virial and 
compressibility pressures; the compressibility pressure, in 
particular, follows the MC results very closely. Surpris- 
ingly, the inclusion of the second-order bridge diagram 
increases the deviation of the pressure from MC at high 
densities. The same figure shows that the m = 2 Kerr 
coefficient agrees more closely with the simulation results 
if bridge corrections are included, but the change is less 
clear for m > 2. 



B. Cavity correlation function 



of the two particles, 1 and 2, and thus may be expanded 
in terms of Legendre functions of U\ ■ u 2 . Using the 
spherical harmonic addition theorem and comparing the 
results with Eq. I|18l) , one finds that y m nx (0) is zer0 unless 
m = n. Furthermore £/ mmx (0) = {-) x y mm o(ty- It may 
be seen from the figures that our calculated functions 
obey this condition to within statistical error. These 
properties result from the fact that the cavity function is 
well-behaved at r = and similar conditions exist for the 
components of the direct correlation function and bridge 
function at r = 0. 

Secondly we note that it has been shown for hard 
spheres that the cavity function at r = is related to 
the excess chemical potential of the fluid, whilst the gra- 
dient of the cavity function at r = is related to the 
pressure [3^, |3|| . These calculations may be generalised 
for anisotropic hard bodies and we obtain the exact re- 
sults (for axially symmetric particles) 



j/ (it, it, r = 0) = exp(/3/x cx ) 
dy(u, u, r) 



dr 



= --j-P2/(«,«,0) 

r=0 47T 



du\ u\ 



ui «>0 



J du 2 rl(ui,u 2 )g(r C: Ui, u 2 ) 



(32) 



(33) 



where the integral over tii is restricted to the positive re- 
gion of its z component; r c (iti, u 2 ) is the contact distance 
of the two ellipsoids and g{r c , u\, u 2 ) is the contact value 
of the pair distribution function at the given orientation. 
In the special case of hard spheres (i.e. r c = constant), 
Eq. H33J) gives the aforementioned relationship with the 
pressure, but in general, so far as we can see, there is no 
simple connection between the r.h.s. of Eq. (|33|l and any 
thermodynamic property of the fluid. 

We now turn to our numerical results. Selected spheri- 
cal harmonics components of y(l, 2) are shown in Figs. El 
0] for two densities p* = 0.1, 0.3. The most obvious 
conclusion is that both the PY and HNC predictions 
differ greatly from the simulation results as the den- 
sity increases. In general the PY predictions are far too 
small in magnitude, whereas the HNC results are far too 
big. This is particularly evident for the isotropic com- 
ponent, 2/000(7"), at low values of r, where yooo{f) rises 
dramatically. At higher densities HNC and simulation 
coefficients differ by several orders of magnitude (from 
simulation yooo(0) = 2050.1, while from HNC theory 
yooo (0) = 3033586.2 for p* = 0.50). 

The inclusion of the bridge diagrams in the HNC clo- 
sure improves significantly the agreement with MC at 
p = 0.1, see Fig. |31 but theory is still far from simulation 
for p = 0.3, see Fig. 01 



Before presenting our numerical results, it is worth 
considering some exact, analytical properties of the cav- 
ity function at r — 0. Firstly we note that at r = 0, the 
cavity function only depends on the relative orientations 



C. Bridge function 

Shown in Figs. I5l7l are selected bridge function compo- 
nents calculated from simulation, PY and virial expan- 
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FIG. 5: Spherical harmonics components of b(l,2) for p* = 
0.10 found from simulation, PY theory , and the virial expan- 
sion truncated at 2nd (62) and 3rd order (63). 




FIG. 6: p* = 0.2 with same functions as in Fig. [S] 



sion truncated at the second order (62(1,2)) and third 
order (63(1,3)). 

As can be seen the PY 6000 ( r ) is always larger than 
the simulation 6000 (f) by approximately a factor of two. 
This seems to be independent of density. The shape of 
this component, both from simulation and PY theorVjis 
similar to that of 6(r) calculated for simple fluids [a,L|- 
The slope of b 00 o(r) goes toward as r goes to 0. This is 
similar to the behaviour seen for 6(r) for Lennard-Jones 
and soft sphere systems 0, Q, while for the HS fluid 
6(r) approaches r = almost linearly |4^. PY theory 
similarly overestimates the angular b rnnx (r) coefficients. 

Shown in Figs. 15171 are also the bridge function cal- 
culated from simulation and from the first and second 
terms in the diagrammatic expansion. As can be seen 
at the lowest density studied (p* — 0.10) the first-order 
expansion gives reasonable agreement with the simulated 
bridge function components, although they are underesti- 
mated relative to simulation. Adding the second term in 
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FIG. 7: p* — 0.5 with same functions as in Fig. 



the expansion improves the agreement quite considerably. 
At a higher density, p* = 0.20, the agreement is less good, 
with the first-order expansion seriously underestimating 
the coefficients. Again the agreement improves with the 
addition of the second term, in fact this diagrammatic ap- 
proximation of the bridge function is better overall than 
the approximation obtained from the PY equation. As 
the density increases the agreement between MC results 
and the truncated virial expansion worsens, see Fig. \7\ 



VI. CONCLUSIONS 

In this paper we have presented the calculation of 
the pair correlation functions h(l, 2), c(l, 2), y(l, 2), and 
6(1,2) for the spheroid fluid from both simulation and 
IET. The total and direct correlation functions have been 
calculated using methods previously described 0, |2(| . 
The cavity function was calculated from simulation us- 
ing a combination of a direct simulation method and a 
test-particle approach. In order to improve the sampling 
of y(l,2) in the direct simulation approach an umbrella 
sampling scheme using a weight function determined it- 
crativcly during the simulation itself is employed. From 
IET the cavity function is determined directly using the 
approximate closure relations. 

Comparison between simulation and integral equation 
show, as reported before El , reasonable agree- 
ment between the coefficients of the total and direct cor- 
relation functions. However theory predicts the simu- 
lated cavity function poorly, with PY theory underes- 
timating and HNC theory overestimating y(l,2) within 
the overlap region. This error rapidly increases with den- 
sity, leading to, at the highest densities studied, errors of 
several orders of magnitude. 

The bridge function calculated from the truncated 
virial expansion is in good agreement with MC results 
at low density but significant differences appear as the 
density increases. The bridge function calculated from 
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PY theory follows the general shape of the MC results 
but the quantitative agreement is poor. 

To the best of our knowledge this work presents the 
first calculation of both the full bridge and cavity func- 
tions for molecular fluids from simulation. As the ap- 
proximate closure relations used in integral equation the- 
ory correspond to approximations to the bridge function, 
knowledge of its exact form will, hopefully, be of great 
benefit in developing improved theories of molecular flu- 
ids. 
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APPENDIX: WANG-LANDAU SAMPLING 

Consider a system with a property X. The probability 
of finding the system with a particular X = X\ is given 
by a probability distribution p(X). In many cases this 
distribution is peaked around certain values of X, mean- 
ing that in a standard simulation values away from these 
are likely to be poorly sampled. When it is desirable to 
get information about these unlikely states it is common 
to apply a weight function, g(X) — exp(— (3W(X)), that 
changes the standard Metropolis acceptance criteria to 

P (X 1 - X 2 ) = CX P h W*a) - E(X 1 ))} . 

9{X2) 

(A.l) 

The simulated probability distribution p s i m (X) then be- 
comes 

p sim (X)=p(X)g(X). (A.2) 

Ideally the effect of the weight function is to make the 
simulation probability distribution flat, i.e. p s i m (X) = 1, 
which implies 

W(X) = hogp(X), (A.3) 
P 

Of course the W(X) needed to achieve this perfectly flat 
histogram is not known in advance, otherwise the proba- 
bility distribution would also be known in advance, thus 
rendering the actual act of performing the simulation 
somewhat redundant. The problem has then become one 
of determining the weight function needed to produce a 
flat histogram. 

At the start of the simulation the weight function is 
initially set to be constant, i.e. g(X) = 1, W(X) = 0. 



After each attempted MC move X% —> X 2 (made using 
the modified criteria Eq. lA.lfl the weight function for the 
resulting state X 1 / 2 (either Xi or X 2 ) is multiplied by a 
modification factor, 

g(x 1/2 ) - fg(x 1/2 ) 

W(X 1/2 )^W(X 1/2 ) + logf . 

Simultaneously the probability histogram p s i m (Xi/ 2 ) is 
also incremented. This continues with the weight func- 
tion and probability being updated after every attempted 
change in X until probability histogram is flat. The flat- 
ness condition may be defined in several ways and will 
be discussed momentarily. Once this condition has been 
reached the probability histogram is reset to zero and / 
is modified. Typically 

/- V 7 / 
log / -> ^ lo g / ■ 

This then continues until the modification factor becomes 
close to 1 (log / gets close to machine precision) . The 
final p(X) may then be determined from 

p(X)=p siul {X)/g(X)=p siul {X)eiq ? {0W(X)) . (A.4) 

A few general notes on the method are due. First up- 
dating the weight function during the simulation may be 
seen to violate the principle of detailed balance. However, 
this is most severe at the beginning of the simulation. 
As / tends toward 1, the changes in the weight func- 
tion become increasingly small. It has been shown that 
a viable MC scheme need only asymptotically obey de- 
tailed balance [5(j . Additionally once a sufficiently good 
weight function has been determined, the simulation may 
be continued without updating the weight function and 
statistics may be gathered from this j4lj. Secondly as a 
perfectly flat histogram is unlikely to be reached during 
a finite simulation the flatness condition may be seen to 
be somewhat arbitrary. In the first implementations the 
histogram was declared flat when the smallest p s - lnl (X) 
was within a given percentage of the average. However 
it is not impossible to imagine pathological distribution 
(e.g. with a few large narrow peaks) that are far from flat 
but still fulfil this criteria. An alternative is to update 
/ whenever every bin has been visited a minimum num- 
ber of times. While this may appear less rigorous then 
the first method, as the simulation progresses and W(X) 
becomes closer to (1/(3) logp(X) thenp s ; m (A) should be- 
come flat. Additionally this ensures that p s i m (X) has a 
chance to adjust to the new / and avoids any spurious 
early updates. One final point is that the Wang-Landau 
method was originally formulated for systems with dis- 
crete degrees of freedom (specifically the Ising model). 
When X is continuous the probability histogram and 
weight functions are calculated for bins of finite width 
X, X + SX and bin width may become a perturbing 
factor in the results. 
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In the present problem we are interested in the prob- 
ability distribution of a pair of non-interacting parti- 
cles. The variable of interest is the radial separation of 
these particles ri2, which is discretized into bins of width 
Sr = 0.01. As mentioned before the r\i range is divided 
into a set of overlapping windows. A Wang-Landau sim- 



ulation is used to determine the weight function to pro- 
duce a constant p(r 12) within each window. / is updated 
whenever p(r 12) fulfils two criteria: i) the largest differ- 
ence between any bin and the average is less than 10% 
and ii) the smallest values of any bin is 100. 
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